import functools
import operator
import numpy as np
def spiht_bitget(number, pos):
return (number >> pos) & 1
def spiht_descendant(i, j, type, m):
s = m.shape[0]
S = []
index = 0
a = 0
b = 0
while 2*i < s and 2*j < s :
a = i
b = j
mind = np.arange(2 * (a + 1) - 2, 2 * (a + np.exp2(index).astype(int)))
nind = np.arange(2 * (b + 1) - 2, 2 * (b + np.exp2(index).astype(int)))
# Check ranges
mind = mind[mind<s]
nind = nind[nind<s]
S.append(m[np.ix_(mind,nind)].flatten())
index = index + 1
i = 2 * a
j = 2 * b
S2 = functools.reduce(operator.iconcat, S, [])
if type == 1:
S2 = np.delete(S2,np.arange(0,4),axis=0)
r = np.abs(np.array(S2))
r = np.max(r)
return r
def spiht_encode(m, max_bits, level):
block_size = m.size
# Initialization
bitctr = 0
out = 2 * np.ones(max_bits - 14,int)
n_max = np.floor(np.log2(np.abs(np.max(m))))
Bits_Header = 0
Bits_LSP = 0
Bits_LIP = 0
Bits_LIS = 0
# output bit stream header
# image size, number of bit plane, wavelet decomposition level should be
# written as bit stream header.
out[0:3] = [m.shape[0], n_max, level]
bitctr = bitctr + 24
index = 3
Bits_Header = Bits_Header + 24
# Initialize LIP, LSP, LIS
temp = []
bandsize = np.exp2(np.log2(m.shape[0]) - level + 1).astype(int)
temp1 = np.arange(bandsize)
temp = np.tile(temp1, (bandsize, 1)).transpose()
LIP = np.zeros((temp.size, 2), int)
LIP[:, 0] = temp.flatten()
LIP[:, 1] = temp.transpose().flatten()
LIS = np.zeros((temp.size, 3), int)
LIS[:, 0] = LIP[:, 0]
LIS[:, 1] = LIP[:, 1]
LIS[:, 2] = np.zeros(len(LIP[:, 0]))
LSP = []
pstart = 0
pend = bandsize // 2
for i in range(bandsize // 2):
LIS = np.delete(LIS, np.array(range(pstart, pend)), axis=0)
pdel = pend - pstart;
pstart = pstart + bandsize - pdel
pend = pend + bandsize - pdel
n = n_max;
# coding
while bitctr < max_bits:
# Sorting Pass
LIPtemp = LIP
temp = 0
for i in range(LIPtemp.shape[0]):
if (bitctr + 1) >= max_bits:
if (bitctr < max_bits):
out = np.delete(out, -1)
return out
if np.abs(m[LIPtemp[i, 0], LIPtemp[i, 1]]) >= np.exp2(n): # 1: positive; 0: negative
out[index] = 1
bitctr = bitctr + 1
index = index + 1
Bits_LIP = Bits_LIP + 1
sgn = m[LIPtemp[i, 0], LIPtemp[i, 1]] >= 0
out[index] = sgn
bitctr = bitctr + 1
index = index + 1
Bits_LIP = Bits_LIP + 1
LSP.append(LIPtemp[i, :].tolist())
LIP = np.delete(LIP, temp, axis=0)
temp = temp - 1
else:
out[index] = 0
bitctr = bitctr + 1
index = index + 1
Bits_LIP = Bits_LIP + 1
temp = temp + 1
LIStemp = LIS
temp = 0
i = 0
while i < LIStemp.shape[0]:
if LIStemp[i,2] == 0:
if bitctr >= max_bits:
return out
max_d = spiht_descendant(LIStemp[i,0], LIStemp[i,1], LIStemp[i,2], m)
if max_d >= np.exp2(n):
out[index] = 1
bitctr = bitctr + 1
index = index + 1
Bits_LIS = Bits_LIS + 1
x = LIStemp[i,0]
y = LIStemp[i,1]
if (bitctr + 1) >= max_bits:
if bitctr < max_bits:
out = np.delete(out, -1)
return out
if np.abs(m[2*x,2*y]) >= np.exp2(n):
LSP.append([2*x,2*y])
out[index] = 1
bitctr = bitctr + 1
index = index + 1
Bits_LIS = Bits_LIS + 1
sgn = m[2*x,2*y] >= 0
out[index] = sgn
bitctr = bitctr + 1
index = index + 1
Bits_LIS = Bits_LIS + 1
else:
out[index] = 0
bitctr = bitctr + 1
index = index + 1
Bits_LIS = Bits_LIS + 1
LIP = np.vstack((LIP,[2*x,2*y]))
if (bitctr + 1) >= max_bits:
if bitctr < max_bits :
out = np.delete(out, -1)
return out
if np.abs(m[2*x,2*y + 1]) >= np.exp2(n):
LSP.append([2*x, 2*y + 1])
out[index] = 1
bitctr = bitctr + 1
index = index + 1
Bits_LIS = Bits_LIS + 1
sgn = m[2*x,2*y + 1] >= 0
out[index] = sgn
bitctr = bitctr + 1
index = index + 1
Bits_LIS = Bits_LIS + 1
else:
out[index] = 0
bitctr = bitctr + 1
index = index + 1
Bits_LIS = Bits_LIS + 1
LIP = np.vstack((LIP,[2*x,2*y + 1]))
if (bitctr + 1) >= max_bits:
if bitctr < max_bits:
out = np.delete(out, -1)
return out
if np.abs(m[2*x + 1, 2*y]) >= np.exp2(n):
LSP.append([2*x + 1, 2*y])
out[index] = 1
bitctr = bitctr + 1
index = index + 1
Bits_LIS = Bits_LIS + 1
sgn = m[2*x + 1, 2*y] >= 0
out[index] = sgn
bitctr = bitctr + 1
index = index + 1
Bits_LIS = Bits_LIS + 1
else:
out[index] = 0
bitctr = bitctr + 1
index = index + 1
Bits_LIS = Bits_LIS + 1
LIP = np.vstack((LIP, [2*x + 1, 2*y]))
if (bitctr + 1) >= max_bits:
if bitctr < max_bits:
out = np.delete(out, -1)
return out
if np.abs(m[2*x + 1, 2*y + 1]) >= np.exp2(n):
LSP.append([2*x + 1, 2*y + 1])
out[index] = 1
bitctr = bitctr + 1
index = index + 1
Bits_LIS = Bits_LIS + 1
sgn = m[2*x + 1,2*y + 1] >= 0
out[index]= sgn
bitctr = bitctr + 1
index = index + 1
Bits_LIS = Bits_LIS + 1
else:
out[index] = 0
bitctr = bitctr + 1
index = index + 1
Bits_LIS = Bits_LIS + 1
LIP = np.vstack((LIP, [2*x + 1, 2*y + 1]))
if (4*x+2 < m.shape[0]) and (4*y+2 < m.shape[1]):
LIS = np.vstack((LIS, [LIStemp[i,0],LIStemp[i,1],1]))
LIStemp = np.vstack((LIStemp, [LIStemp[i, 0], LIStemp[i, 1], 1]))
LIS = np.delete(LIS, temp, axis=0)
temp = temp - 1
else:
out[index] = 0
bitctr = bitctr + 1
index = index + 1
Bits_LIS = Bits_LIS + 1
else:
if bitctr >= max_bits:
return out
max_d = spiht_descendant(LIStemp[i,0], LIStemp[i,1], LIStemp[i,2], m)
if max_d >= np.exp2(n):
out[index] = 1
bitctr = bitctr + 1
index = index + 1
x = LIStemp[i,0]
y = LIStemp[i,1]
LIS = np.vstack((LIS,[2*x , 2*y , 0]))
LIS = np.vstack((LIS,[2*x , 2*y + 1, 0]))
LIS = np.vstack((LIS,[2*x + 1, 2*y , 0]))
LIS = np.vstack((LIS,[2*x + 1, 2*y + 1, 0]))
LIStemp = np.vstack((LIStemp,[2*x , 2*y , 0]))
LIStemp = np.vstack((LIStemp,[2*x , 2*y + 1, 0]))
LIStemp = np.vstack((LIStemp,[2*x + 1, 2*y , 0]))
LIStemp = np.vstack((LIStemp,[2*x + 1, 2*y + 1,0]))
LIS = np.delete(LIS, temp, axis=0)
temp = temp - 1
else:
out[index] = 0
bitctr = bitctr + 1
index = index + 1
Bits_LIS = Bits_LIS + 1
i = i + 1
temp = temp + 1
# Refinement Pass
temp = 0
value = np.floor(np.abs(np.exp2(n_max - n + 1) * m[LSP[temp][0],LSP[temp][1]])).astype(int)
while (value >= np.exp2(n_max + 2) and (temp < len(LSP))):
if bitctr >= max_bits:
return out
s = spiht_bitget(value, int(n_max + 1))
out[index] = s
bitctr = bitctr + 1
index = index + 1
Bits_LSP = Bits_LSP + 1
temp = temp + 1
if temp < len(LSP):
value = np.floor(np.abs(np.exp2(n_max - n + 1) * m[LSP[temp][0],LSP[temp][1]])).astype(int)
n = n - 1
return out
def spiht_decode(bits):
# Initialization
m = np.zeros((bits[0],bits[0]),np.double)
n_max = bits[1]
level = bits[2]
ctr = 3
# Initialize LIP, LSP, LIS
temp = []
bandsize = np.exp2(np.log2(bits[0]) - level + 1).astype(int)
temp1 = np.arange(bandsize)
temp = np.tile(temp1, (bandsize, 1)).transpose()
LIP = np.zeros((temp.size, 2), int)
LIP[:, 0] = temp.flatten()
LIP[:, 1] = temp.transpose().flatten()
LIS = np.zeros((temp.size, 3), int)
LIS[:, 0] = LIP[:, 0]
LIS[:, 1] = LIP[:, 1]
LIS[:, 2] = np.zeros(len(LIP[:, 0]))
LSP = []
pstart = 0
pend = bandsize // 2
for i in range(bandsize // 2):
LIS = np.delete(LIS, np.array(range(pstart, pend)), axis=0)
pdel = pend - pstart;
pstart = pstart + bandsize - pdel
pend = pend + bandsize - pdel
# coding
n = n_max
while ctr < bits.shape[0]:
# Sorting Pass
LIPtemp = LIP
temp = 0
for i in range(LIPtemp.shape[0]):
if ctr >= bits.shape[0] - 1: # addtl cond., prevents counter overflow
return m
if bits[ctr] == 1:
ctr = ctr + 1
if bits[ctr] > 0:
m[LIPtemp[i,0], LIPtemp[i,1]] = np.exp2(n) + np.exp2(n - 1)
else:
m[LIPtemp[i,0], LIPtemp[i,1]] = -1*np.exp2(n) - np.exp2(n - 1)
LSP.append(LIPtemp[i, :].tolist())
LIP = np.delete(LIP, temp, axis=0)
temp = temp - 1
ctr = ctr + 1
temp = temp + 1
LIStemp = LIS
temp = 0
i = 0
while i < LIStemp.shape[0]:
if ctr >= bits.shape[0]:
return m
if LIStemp[i,2] == 0:
if bits[ctr] == 1:
ctr = ctr + 1
x = LIStemp[i,0]
y = LIStemp[i,1]
if ctr >= bits.shape[0] - 1: # addtl cond., prevents counter overflow
return m
if bits[ctr] == 1:
LSP.append([2*x,2*y])
ctr = ctr + 1
if bits[ctr] == 1:
m[2*x, 2*y] = np.exp2(n) + np.exp2(n - 1)
else:
m[2*x, 2*y] = -1*np.exp2(n) - np.exp2(n - 1);
ctr = ctr + 1
else:
LIP = np.vstack((LIP,[2*x,2*y]))
ctr = ctr + 1
if ctr >= bits.shape[0] - 1: # addtl cond., prevents counter overflow
return m
if bits[ctr] == 1:
ctr = ctr + 1
LSP.append([2*x, 2*y + 1])
if bits[ctr] == 1:
m[2*x, 2*y+1] = np.exp2(n) + np.exp2(n - 1)
else:
m[2*x, 2*y+1] = -1*np.exp2(n) - np.exp2(n - 1)
ctr = ctr + 1
else:
LIP = np.vstack((LIP,[2*x, 2*y + 1]))
ctr = ctr + 1
if ctr >= bits.shape[0] - 1: # addtl cond., prevents counter overflow
return m
if bits[ctr] == 1:
ctr = ctr + 1
LSP.append([2*x + 1, 2*y])
if bits[ctr] == 1:
m[2*x + 1, 2*y] = np.exp2(n) + np.exp2(n - 1)
else:
m[2*x + 1, 2*y] = -1*np.exp2(n) - np.exp2(n - 1)
ctr = ctr + 1
else:
LIP = np.vstack((LIP,[2*x + 1, 2*y]))
ctr = ctr + 1
if ctr >= bits.shape[0] - 1: # addtl cond., prevents counter overflow
return m
if bits[ctr] == 1:
ctr = ctr + 1
LSP.append([2*x + 1, 2*y + 1])
if bits[ctr] == 1:
m[2*x + 1, 2*y + 1] = np.exp2(n) + np.exp2(n-1)
else:
m[2*x + 1, 2*y + 1] = -1*np.exp2(n) - np.exp2(n-1)
ctr = ctr + 1
else:
LIP = np.vstack((LIP,[2*x + 1, 2*y + 1]))
ctr = ctr + 1
if (4*x+2 < m.shape[0]) and (4*y+2 < m.shape[1]):
LIS = np.vstack((LIS, [LIStemp[i,0],LIStemp[i,1],1]))
LIStemp = np.vstack((LIStemp, [LIStemp[i, 0], LIStemp[i, 1], 1]))
LIS = np.delete(LIS, temp, axis=0)
temp = temp - 1
else:
ctr = ctr + 1
else:
if bits[ctr] == 1:
x = LIStemp[i,0]
y = LIStemp[i,1]
LIS = np.vstack((LIS,[2*x , 2*y , 0]))
LIS = np.vstack((LIS,[2*x , 2*y + 1, 0]))
LIS = np.vstack((LIS,[2*x + 1, 2*y , 0]))
LIS = np.vstack((LIS,[2*x + 1, 2*y + 1, 0]))
LIStemp = np.vstack((LIStemp,[2*x , 2*y , 0]))
LIStemp = np.vstack((LIStemp,[2*x , 2*y + 1, 0]))
LIStemp = np.vstack((LIStemp,[2*x + 1, 2*y , 0]))
LIStemp = np.vstack((LIStemp,[2*x + 1, 2*y + 1,0]))
LIS = np.delete(LIS, temp, axis=0)
temp = temp - 1
ctr = ctr + 1
i = i + 1
temp = temp + 1
# Refinement Pass
temp = 0
value = m[LSP[temp][0], LSP[temp][1]]
while (np.abs(value) >= np.exp2(n + 1) and (temp < len(LSP))):
if ctr >= bits.shape[0] - 1: # addtl cond., prevents counter overflow
return m
value = value + np.float_power(-1, bits[ctr] + 1) * np.float_power(2,n - 1) * np.sign(m[LSP[temp][0], LSP[temp][1]])
m[LSP[temp][0], LSP[temp][1]] = value
ctr = ctr + 1
temp = temp + 1
if temp < len(LSP):
value = m[LSP[temp][0], LSP[temp][1]]
n = n - 1
return m
